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Time-Dependent Density Functional Theory (TDDFT) has recently been extended to describe many-body 
open quantum systems (OQS) evolving under non-unitary dynamics according to a quantum master equation. 
In the master equation approach, electronic excitation spectra are broadened and shifted due to relaxation and 
dephasing of the electronic degrees of freedom by the surrounding environment. In this paper, we develop a 
formulation of TDDFT linear-response theory (LR-TDDFT) for many-body electronic systems evolving under 
a master equation, yielding broadened excitation spectra. This is done by mapping an interacting open quantum 
system onto a non-interacting open Kohn-Sham system yielding the correct non-equilibrium density evolution. 
A pseudo-eigenvalue equation analogous to the Casida equations of usual LR-TDDFT is derived for the Redfield 
master equation, yielding complex energies and Lamb shifts. As a simple demonstration, we calculate the 
spectrum of a atom in an optical resonator interacting with a bath of photons. The performance of an 
adiabatic exchange-correlation kernel is analyzed and a first-order frequency-dependent correction to the bare 
Kohn-Sham linewidth based on Gorling-Levy perturbation theory is calculated. 

PACS numbers: 



I. INTRODUCTION 

Due to its attractive balance between accuracy and effi- 
ciency, time-dependent density functional theory (TDDFT) 
has seen a tremendous growth of applications in recent years. 
These range from optical properties of molecules, clusters and 
solids, to optimal control theory and real-time dynamics of 
species in intense laser fields L1-;6J- TDDFT has been partic- 
ularly successful at calculating optical response properties of 
electronic systems in the linear response regime [7]. In most 
quantum chemical codes, excitation energies and oscillator 
strengths are extracted by solving a pseudo-eigenvalue equa- 
tion, originally formulated by Casida fs*]. The Casida equa- 
tions are derived by considering the linear density response 
of an interacting system and corresponding non-interacting 
Kohn-Sham system, both undergoing unitary evolution. If the 
Casida equations are solved using the ubiquitous adiabatic ap- 
proximation (ATDDFT) within a discrete basis set, the result- 
ing eigenvalues are real. This gives rise to a discrete absorp- 
tion spectrum of delta function peaks. 

In experimentally observed spectra, line broadening arises 
from a variety of different mechanisms, several of which 
have been explored akeady within LR-TDDFT. In extended 
systems, relaxation and dephasing due to electron-electron 
scattering is well captured us ing no n-adiabatic and cuiTent- 
density dependent functionals H-IJ]. In finite systems, decay 
of resonant states due to coupling with the continuum gives 
rise to finite linewidths. The ability of DFT and TDDFT to 
capture lineshape parameters and widths of resonances has 
been discussed in 



Another important broadening mechanism arises from re- 
laxation and dephasing of electronic degrees of freedom by 
a classical or bosonic bath such as photons, phonons or im- 
purities. For extended systems, this situation was consid- 
ered in 1 3- 18 1, where linewidths of intersubband plasmon 



excitations were well captured by combining the Vignale- 
Kohn functional |9] to account for electron-electron scattering 
with the memory function formalism for electron-phonon and 
electron-impurity scattering. For atomic and molecular sys- 
tems, the theory of open quantuni systems within the master 
equation approach is often used Several important 

examples include vibrational relaxation of molecules in liq- 
uids and solid impurities fls - 27 ] , cavity quantum electrody- 
namics (QED) 128-37], photo-absorption of chromophores in 
a protein bath [32-35], single-molecule transport [36-40] and 
exciton transport |l4li-(44t|. In all of these examples, even with 
simple system-bath models, the exact solution of the master 
equation for the reduced dynamics of the many-body elec- 
tronic system is computationally intractable. Therefore, open 
quantum systems TDDFT (OQS -TDDFT) offers an attractive 
approach to the many-body open-systems problem. 

Several important steps toward the formulation of OQS- 
TDDFT have recently been made, with the focus on real- 
time dynamics. In [36], a Runge-Gross theorem was es- 
tablished for Markovian master equations of the Lindblad 
form. A scheme in which the many-body master equation 
is mapped onto a non-interacting Kohn-Sham master equa- 
tion was proposed for application to single-molecule trans- 
port. In 145114611 . the Runge-Gross theorem was extended to 
arbitrary non-Markovian master equations and a Van Leeuwen 
construction was established, thereby proving the existence of 
an OQS-TDDFT Kohn-Sham scheme iH. In H, it was 
shown that the original interacting open system dynamics can 
be mapped onto either a non-interacting open Kohn-Sham sys- 
tem, or a non-interacting closed (unitarily evolving) Kohn- 
Sham system. A different formulation of OQS-TDDFT based 
on the stochastic Schrodinger equation has also been devel- 
oped EsUlS]. 

The goal of the present manuscript is to formulate the linear 
response version of (OQS-TDDFT) within the master equa- 
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tion approach. This provides a framework in which environ- 
mentally broadened spectra of many-body electronic systems 
can be accessed in an ab initio way using TDDFT, especially 
when combined with microscopically derived master equa- 
tions. We use the scheme discussed in ll36l, |45], in which 
the interacting OQS can be mapped onto a non-interacting 
open Kohn-Sham system, yielding the same density response. 
This scheme is better suited to response theory than the closed 
Kohn-Sham scheme also discussed in |45], since relaxation 
and dephasing is already accounted for in the Kohn-Sham 
system. The unknown (OQS-TDDFT) exchange-correlation 
functional only needs to correct the relaxation and dephasing 
in the Kohn-Sham system to that of the interacting system, 
rather than needing to explicitly account for the entire effect 
of the environment. However, the closed Kohn-Sham scheme 
is better suited for real-time dynamics, since one only needs 
to propagate a set of equations for the Kohn-Sham orbitals as 
in usual TDDFT. This is in contrast to the open Kohn-Sham 
scheme, in which A^^ — 1 equations are propagated for the ele- 
ments of the density matrix, with N being the dimensionality 
of the Hilbert space. 



II. GENERAL FORMULATION OF OQS-TDDFT LINEAR 
RESPONSE THEORY 

A. Linear response of interacting open quantum systems 

Our formulation of the interacting OQS density-density re- 
sponse function parallels that used in 15211 for calculating spin 
susceptibilities (see also [53', The starting point is the 

unitary evolution for the full density matrix of the system and 
the reservoir (we use the terms "reservoir" and "bath" inter- 
changeably throughout) , 

at I 

where ^(f ) is the Liouvillian superoperator for the full sys- 
tem and reservoir dynamics. The full Hamiltonian is given by 

H{t)=Hs{t)+HR + V. (2) 

Here, 



The paper is organized as follows. In Section II, we formu- 
late the most general OQS-TDDFT linear response equations 
for arbitrary non-Markovian master equations with initial cor- 
relations. In section III, we make the treatment more specific 
by focusing on the Redfield master equation. We also derive 
Casida-type equations whose solution yields the environmen- 
tally broadened absorption spectrum. The solutions to these 
equations are complex, with the real part of the frequency 
yielding the location of absorption peaks and the imaginary 
part yielding the linewidths. In Section IV, we apply the for- 
malism developed in Section III to a atom in an opti- 
cal resonator setup evolving under the Redfield master equa- 
tion. Section V analyzes the performance of using an adia- 
batic functional (OQS-ATDDFT) in solving the OQS-TDDFT 
Casida equations derived in section III. To a large degree, 
OQS-ATDDFT is seen to provide a reliable correction to the 
location of absorption peaks while leaving the linewidths un- 
changed. A frequency-dependent functional yielding a correc- 
tion to the OQS-ATDDFT linewidth based on Gorling-Levy 
(GL) perturbation theory is then calculated and analyzed. In 
section VI a discussion and outlook is provided. 



m) = L Vr + £ + Lvv.« (r„ f ), (3) 

is the Hamiltonian of the electronic system of interest in an 
external potential V(.a(r,f). This potential generally consists 
of a static external potential due to the nuclei and an external 
driving field coupled to the system such as a laser field. The 
system-bath coupling, V, is at this point arbitrary, but we will 
discuss specific forms later. V acts in the combined Hilbert 
space of the system and reservoir and so it couples the two 
subsystems. Hr is the Hamiltonian of the reservoir, assumed 
to have a dense spectrum of eigenstates. The density of states 
of Hr determines the structure of reservoir correlation func- 
tions, whose time-scale in turn determines the reduced system 
dynamics. 

Defining the reduced density operator for the electronic sys- 
tem alone by tracing over the reservoir degrees of freedom, 

Psit)^TrR{p{t)}, (4) 
one arrives at the formally exact quantum master equation. 



Atomic units in which e — h — mg — I are used through- 
out. This also implies that the speed of light in vacuum is 
given by c = 137. For generality, we have formulated most of 
the theory by considering linear response from an equilibrium 
state at finite temperature. For atoms and molecules, it will 
generally be sufficient to take the zero temperature limit and 
consider linear response from the ground-state. 



j^psit) - 'i[Hsit),psit)]+ J' dTZ{t~T)ps{T)+'i'{t). 

(5) 

Here, E(f — t) is the memory kernel and ^{t) arises from 
initial correlations between the system and its environment. 
It is referred to as the inhomogeneous term. The above equa- 
tion is still formally exact, as ps{t) gives the exact expectation 
value of any observable depending only on the electronic de- 
grees of freedom. In practice, however, approximations to S 
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and are required. Of particular importance in TDDFT is the 
time-dependent electronic density, 



n(r,t) =e'"'n{r)e 



„-tHt 



(15) 



«(r,f)-rrs{ps(f)«(r)}, 



(6) 



where «(r) = L/^5(r — f,). We now assume that for t < to, 
the external potential is time-independent while for t > to a 
weak perturbing field is applied, i.e. 



t <to,Vext{r,t) = vw(r) (7) 
t >fo,vw(r,f) = vw(r) + 5vvvr(r,f)- (8) 

For f < fo, the entire system and environment is in thermal 
equilibrium described by the canonical density operator 



is the operator generating the electronic charge density in 
the Heisenberg picture with respect to the full Hamiltonian for 
t < fo- Rearranging terms under the trace operation in Eq.[T4l 
the density-density response function can be written as 

X„„ir,r'; CO) = lim i / dte-""'-''Trs{n(r)p"sir\t)}. (16) 

£^+0 Jo 

Here, p'^{r,t) is an operator acting in the electronic system 
Hilbert space, which obeys the same equation of motion as the 
reduced system density operator (Eq.|5]) 



-pH 



(9) 



(17) 

subject to the initial condition 



where /3 = is the inverse temperature. The reduced 
equilibrium density operator of the electronic system is then 
given by 



p"s{r,0) = Tru{[n{r),p'"^]}. 



(18) 



Ps 



eq 



(10) 



Carrying out the Fourier transform in Eq. [161 one arrives 
at the formally exact expression for the open-systems density- 
density response function in Liouville space. 



In Eq. IHand Eq. [TOl // = + + V is the full Hamilto- %„„ (r, r'; to) = iTrs { «(r) 



1 



nian for r < fo and = - i I^i V2 + y!^^. ^ + J]. v,.« (r,) 
is the static Hamiltonian of the electrons in the absence of the 
external perturbation. 

For t > to, the perturbing field is switched on and the sys- 
tem density operator subsequently evolves under the master 
equation given in Eq. |5] The electronic density evolution to 
first-order in the perturbing field is then given by 



t <fo,n(r,f)-'/nr) (H) 
f >fo,n(r,f) =n''''(r) + 5«(r,f). (12) 

Here, rf^{r) = Trs{pl^h{ry} is the equilibrium electron 
density and 



5n(r,f) = j d'r' j dt' x„n{r,t;r' ,t')8veAr\t') (13) 

is the lineal" density response. ;^„„(r,f;r',f') is the density- 
density response function. Its Fourier transform to the fre- 
quency domain is given by 

poo 

Xnn{r,r';(0)^ lim i / dte-""-''Trs+R{[n{r,t),n{r'W'^}, 
e^+O Jo 

(14) 

where 



■(p,"(r',0)+'P(«))j 
(19) 

J^s is the Liouvillian for the system Hamiltonian for t < to, 
defined by it's action on an arbitrary operator O by J/fgO = 
[Hs,0]. It is readily verified that Eq.[T9l reduces to the usual 
expression for the density-density response function of an 
isolated system when Z{co) — 0, ^(o)) = and p^(r,0) ~ 
[«(r),pf]. 

The absorption spectrum can be extracted by taking the 
imaginary part of ;if„„(r,r'; ©) in Eq.[T9] For an isolated sys- 
tem with a discrete spectrum, this is given by a sum over 
weighted delta function peaks. For an open system as in 
Eq.[T9l S(a)) and ^(o)) in principle give rise to the exact com- 
plicated broadened and shifted spectrum, due to relaxation 
and dephasing of the electronic degrees of freedom by the en- 
vironment. In practice, however, even with simple approxima- 
tions to S((b) and ^{co), the exact form of 3m [;^„„(r,r', O))] 
is not exactly known, since it refers to a many-body response 
function. In the next two subsections, we consider an open 
Kohn-Sham system, formally yielding the exact density re- 
sponse of the original interacting open system. In an open- 
systems TDDFT framework, the exact spectrum of Eq. [T9]is 
obtained by correcting the open Kohn-Sham spectrum via an 
exchange-correlation kernel. The kernel must take into ac- 
count not only the electron-electron interaction contained in 
but must also correct the interaction of the system with 
the bath, described by £((») and ^'{co). 
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B. The open Kohn-Sham system 



It was proven in II45I1 . that for a master equation of the form 
given in Eq.|5] there exists a unique, non-interacting and open 
Kohn-Sham system, whose system density operator evolves 
under the master equation 

(20) 

such that the time-dependent density is obtained from 



n{r,t) = Tr{pl\t)n{r)} (21) 

for all times. i/*^"(f) = L^i/if (0- where the Kohn-Sham 
Hamiltonian is given by 



h'-\r,t) 



-vks{r,t). 



(22) 



Here, v^.^ is a local, multiplicative, one-body potential 
which drives the open Kohn-Sham system in such a way that 
the true density of the original interacting open system is re- 
produced for all times. In analogy to usual TDDFT, the Kohn- 
Sham potential is partitioned as 



in Eq.|24]is automatically removed. The equilibrium density, 
rf^{r), is obtained by solving the Kohn-Sham-Mermin equa- 
tions 115 in 



vl?[«](r) 



(/)i(r) = ei0/(r). 



The Kohn-Sham-Mermin potential is partitioned as 



V J [n] (r) = Vext [n] (r) + v/, [«] (r) + 



(25) 



(26) 



5«(r) 

where Fvc[n] is the exchange-correlation contribution to the 
free energy. After solving Eq. |25j the equilibrium Kohn- 
Sham-Mermin density operator is obtained by populating the 
orbitals according to 



1^; = I/i#->(0-l, 
where /, are Fermi-Dirac occupation numbers 

1 



(27) 



(28) 



+ 1 

Denoting (r|0,) = ^,(r), the one-particle Kohn-Sham 
Mermin density matrix is 



vfa(r,f) = v,v,(r,r) + v/,(r,f) +y;r"(r,0, (23) 

where v/,(r,r) is the Hartree potential and the unknown 
functional v"c'^"(r,f) accounts for electron-electron interac- 
tion within the system as well as interaction between the sys- 
tem and bath. In general, 

vr"(r,f) = vT"{r,t)[nX^'','V,^''^\ps{0),p'\Q)]. (24) 

Formally, the open-systems exchange-correlation potential 
is a functional not only of the density, but also of the memory 
kernel, inhomogeneous term and initial state of both the inter- 
acting and Kohn-Sham systems. It has been shown in usual 
TDDFT of closed systems, that initial state dependence can 
be absorbed as dependence on the history of the density and 
vice versa [,5 6;] . Interestingly, in the theory of open quantum 
systems, it is possible to absorb the inhomogeneous term ^, 
into the memory kernel E |55]. This raises the possibility that 
v"?"" may be a functional only of n, E and S'^', but a more rig- 
orous study of this will be done in future work. For notational 
convenience, we suppress the explicit functional dependence 
of v"c on these quantities, although it is implied unless oth- 
erwise stated. 

In general, E*^* and ^'^^ can be chosen to simplify v?c'^"(r,f ) 
as much as possible, although with some restrictions ll45ll and 
consistency conditions between E*^* and 'P*^* fssl]. 

If the system is started in an equilibrium state, as is typically 
the case in linear response theory, the initial state dependence 



{r\t:,\v') ^ 7(r,r') = Y^Mt{vmv'). (29) 

(=1 

The equilibrium density is then obtained by taking the di- 
agonal elements in real space, 

«"'(r) = 7(r,r) = f; /,#,-(r)|2. (30) 

(=1 

For the open Kohn-Sham scheme to be useful practically, 
Vxt^inKr,?), S*^* and should be constructed so that the 
following conditions are satisfied: 

1) E*^* and should not induce correlations between 
non-interacting electrons as the Kohn-Sham system evolves. 
This ensures that the N-body Kohn-Sham density matrix in 
Eq. |20]can be traced over N-1 electron coordinates to arrive 
at a closed equation of motion for the Kohn-Sham reduced 
1 -particle density matrix. Physically, this is expected since 
most reasonable bath models couple to the electronic system 
through one-body operators 1 2 1 , 6 1 ] . 

2) The equation of motion for the non-equilibrium open 
Kohn-Sham reduced 1 -particle density matrix should have the 
Kohn-Sham-Mermin density matrix as its stationary-point so- 
lution. This ensures that the system thermalizes to the cor- 
rect equilibrium density. This also means that at equilibrium, 
Vx^"°"[«](r,f) should reduce to ■ Although this is auto- 
matically satisfied by using v^^ [n] as an adiabatic approxima- 
tion, it might not be satisfied by more sophisticated approxi- 
mations with memory dependence 156145911 . 
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C. Linear response of the open Kohn-Sham system 

Returning to linear response, for f < f o the Kohn-Sham sys- 
tem is in thermal equilibrium with its environment at inverse 
temperature j3, described by Eq. |25] At t — to, the pertur- 
bation 5vv;rt(r,f) is switched on and the Kohn-Sham system 
subsequently evolves according to Eq. |20] The Kohn-Sham 
potential is given by 

f <fo,Vfcv(r,r)=v^![n](r) (31) 
t >fo,Vfa(r,r) = v'^n«](r) + ^^A-.v(r,f), (32) 

where 5vfa(r,f) = dve,,{r,t) + 5vh[n]{r) + 5vT"[n]{r,t). 
Due to Eq. (|2TI ). the exact linear density response of Eq. (fT2] i 
is obtained through 



master equation in the next section. Since correlation between 
the open Kohn-Sham system and reservoir is already partially 
captured through £'^*(a)) and the bai-e Kohn-Sham 

absorption spectrum extracted from 3m [x^*(a),r,r')] is al- 
ready broadened and shifted. The functional fxf™[n'^'^] has 
the task of correcting the spectrum extracted from Xm 1° f^hat 
of the interacting j,,,,, incorporating both the usual electron- 
electron correlation in closed-systems TDDFT as well as addi- 
tional system-bath correlation. In general, the memory kernel 
E(a)) may give rise to a very complicated non-analytic struc- 
ture of Xnii in the lower half of the complex plane. However, 
for Markovian master equations, it will be seen that the pole 
structure of Xnn in the discrete part of the spectrum consists 
of simple poles in the lower complex plane, shifted by a fi- 
nite amount off of the real axis. In such cases, it might be 
reasonable to approximate f"c"^" as 



5«(r,f) = / d'r' / dt'x^,:Xr,fy,t')8vUr',t'). (33) /«P-[„^''](r,r'; «) 



Here, Xn^„{i^,t;r' ,t') is the density-density response function 
of the open Kohn-Sham system. It's Fourier transform to the 
frequency domain is given by 



X^:{co,ry) = iTr\n{r) 



(0 + Ji'k,-iZk'<{(o) 



([«(r'),pf(0)]+^'"(®)) ,(34) 



o«(r)on(r ) 

(37) 

Here, the first term is just the adiabatic contribution to the 
exchange-correlation kernel and the second term is an in gen- 
eral frequency-dependent and imaginary correction. This is 
attractive, since we can take advantage of the usual good per- 
formance of adiabatic TDDFT in describing the location of 
absorption peaks, and attempt to build functionals that go be- 
yond the adiabatic approximation to account for line broaden- 
ing and lamb shifts. This strategy will be discussed further in 
section V. 



where is the Liouvillian for the equilibrium Kohn- 
Sham-Mermin Hamiltonian. Since the system is in the equi- 
librium state at f = 0, p^^XO) must yield the equilibrium den- 
sity, implying that it reduces to the Kohn-Sham-Mermin den- 
sity matrix when traced over N-1 electron coordinates. We 
now define the open-systems exchange-correlation kernel in 
analogy to usual TDDFT for closed-systems by. 



fOpe,u„eqy^ gvT" [«] (r, 0)) , 

// [n^](r,r,a)) = —g--^^;j-^^ (35) 

which is a functional of the equilibrium density. As in 
TDDFT for closed systems, the interacting and Kohn-Sham 
response functions are related through a Dyson-like equation. 



;)^„„(a),r,r') = Z™(a),r,r') 

+ / a'Wz„'^(«,r,y){^^+/jr"K1(y,y';«)} 

X X,m{(o,y'y). (36) 



X,,*, being much simpler then the original interacting Xnn, 
can readily be constructed from the orbitals and eigenvalues in 
Eq.lZSland approximations to E*^^ (a)) and ^'^* (a)) in terms of 
these quantities. This will be done explicitly for the Redfield 



III. LR-TDDFT FOR THE REDFIELD MASTER EQUATION. 

In section II, we formulated LR-TDDFT for a very general 
class of master equations. In this section, we make the dis- 
cussion more specific by invoking the Markov approximation 
and second Born approximation in the system-bath interac- 
tion, to arrive at the microscopically-derived Redfield master 
equation ||32> 61, 78 -81.1 . Since the Redfield equations are rig- 
orously obtained without phenomenological parameters, they 
are amenable to an ab initio theory such as TDDFT. Although 
we focus on Redfield theory here, the generalization of our 
formalism to other Markovian master equations can be made 
with small modifications. Finally, we discuss how it is possi- 
ble to extract the absorption spectrum of a many-body system 
evolving under the Redfield equations directly within OQS- 
TDDFT. This is done by formulating Casida-type equations 
yielding complex eigenvalues due to coupling with the bath. 



A. The Markov approximation and the Redfield master 
equation. 

The Markov approximation describes a situation in which 
the bath correlation functions decay on an infinitely fast time- 
scale relative to the thermalization time of the system L6lll78ll . 
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As a result, the bath has no memory and the memory kernel is For a detailed derivation of the Redfield equations see 11781] . 

time-local In Eq.l44l 



S(r-T)oc J'5(f-T). 



(38) 



Additionally, this implies that the initial density operator is 
a tensor product of a density operator in the system space with 
the equilibrium density operator of the bath 



Sab^N J d^r J d^r2...d^rNW*,i^^^2,--rN)Si-,r)\i/h(r,r2,...rN) 

(45) 

are matrix elements of S{p,r) between system many -body 
wavefunctions and 



p(O)=p,(O)0. 



TrR{e-P"R} 



(39) 



As a result of Eq. [39] the system and environment have no 
initial correlations, and 



l'(f) = 0. (40) 
The master equation (Eq.|5]l then takes the simple form. 



(Oab — Ea — El, 



(46) 



are system many-body excitation energies. g{T) are bath 
correlation functions given by 



where 



(t) = TrR{R{T)R{Q)}, 



(47) 



(48) 



d_ 
dt 



(41) 



If the system Hamiltonian is time-independent, Eq. [41] is 
written in a basis of eigenstates of Hs as: 



^Pab{t) = -l(OahPab{t) + F RabcdPcd{t) 



(42) 



Here, cOab — Ea—Ei, are many -body transition frequencies 
of Hs and R^hcd ^re matrix elements of ^ in this basis. So 
far our discussion applies to any Markovian master equation. 
To obtain the Redfield equations, we further assume that the 
system-bath coupling has a biUnear form 



V^~S(g>R, 



(43) 



where R is an operator in the reservoir Hilbert space which 
couples to a local one-body operator S — [L/=i •^(Pfir,)] in 
the system Hilbert space. This form of the system-bath cou- 
pling is very general and can apply to electron-phonon cou- 
pling in molecules and solid impurities, but also momentum 
dependent couplings which are relevant for instance in laser 
cooling, brownian motion in liquids or dissipative strong field 
dynamics 11551 roll 18211 . The Redfield tensor is then derived 
by performing second-order perturbation theory in V, and is 
given explicitly by 



Rabcd = - / dT 

Jo 



B. Linear response of a many-body system evolving under the 
Redfield master equation. 

Since we consider linear response from the equilibrium 
state, the initial density matrix for the system is given by 



PsiO) 



-PHs 



(49) 



Trs{e-P"s}' 

and the density-density response function in Eq.[T9lreduces 



to 



X„„{r,r';co) = iTrs< «(r) 



1 



CO+J^'s-l. 



^[«(r'),P5(0)] 



(50) 

Inserting a complete set of eigenstates of Hs in Eq. |50l a 
sum over states expression for the density-density response 
function is given by, 

X™(r,r';a)) 

Y rl 0} + CO,b + lRabab 

{a\n{r')\b){b\n{r)\a}\ 



(0 - COab + iR 



(51) 



abab 



dhdY^SanSnce'^^-'-SacSdbe'"^" 



Here, P(£'„ 



are equilibrium occupation probabili- 



dx 



ties of the various many-body states. 

By hermiticity of the density matrix, it can be readily ver- 
^^^-jified that = Rbaba- We can separate the real and imagi- 
nary parts of Rahab as 



7 



abah 



- (lb " 



lA. 



■ab- 



(52) 



From the pole structure of Eq. [ST] we see that Tab corre- 
sponds to an imaginary part of the energy of the transition 
(Oub, giving rise to a finite lifetime, while A^^ is a Lamb shift 
of the real part of the energy. The effect of the Redfield tensor 
is to shift the poles of the density-density response function 
by a finite amount into the lower half of the complex plane. 



C. The Markovian Kohn-Sham-Redfield equations 



dr 



Jo 
Here, 



kl 



SikY.ShnSmje"^"''' - SikSlje" 



= e< - £j 



(57) 



(58) 



are bare Kohn-Sham-Mermin transition frequencies and 



We now discuss the properties of the open Kohn-Sham sys- 
tem for the Redfield master equation. As discussed in section 
II B, there is some freedom in the construction of the Kohn- 
Sham dissipative superoperator and corresponding exchange- 
correlation potential. In this section, we choose a very natural 
form for the Kohn-Sham superoperator, whose construction is 
discussed below. 

We consider an open Kohn-Sham system evolving under a 
Markovian master equation 



dt 



pf{t) = -i[H'^\t),pf{t)] +^/'pf{t), (53) 



'ks ^ ks ( 



which reproduces the exact density evolution of the inter- 
acting Markovian master equation in Eq. [41] To satisfy con- 
dition 1 in section II B, we must have 



J*^(ri,r2,...r^)^X;r^-^-(r,), 



(54) 



i.e. the N-body Kohn-Sham dissipative superoperator is 
a sum of one-body superoperators acting on each coordinate 
separately. This also follows very naturally from the assumed 
one-body nature of the system-bath interaction in Eq. |43] 
Since //*^ (f ) = Yli=\ also a sum of one body terms, we 
can trace both sides of Eq. |53]over N-1 electron coordinates 
and arrive at a closed equation of motion for the Kohn-Sham 
1 -particle reduced density matrix. 



d'x^:{x)S{-,x)^j{x) 



(59) 



are matrix elements of the system-bath coupling operator 
between Kohn-Sham-Mermin orbitals. 

Eq. |56] has a number of desirable properties. First, it's 
stationary point solution is the Kohn-Sham-Mermin density 
matrix, Eq. |27] if h^'^it) reduces to the Kohn-Sham-Mermin 
Hamiltonian when evaluated on the equilibrium density. This 
ensures that condition 2 in section II B is satisfied. It also sat- 
isfies detailed balance as well as most other properties of the 
usual many-body Redfield equations, but in terms of Kohn- 
Sham-Mermin quantities. Also, the tensor r^jj.^ has a simple 
form and can be constructed explicitly in terms of orbitals and 
eigenvalues obtained in an equilibrium-state Kohn-Sham cal- 
culation. The potential v.te "'(f) contained in /!*^* (f) will in gen- 
eral be a functional of y^^ and ^ as well as the time evolving 
density. 



D. Linear response of the Kohn-Sham-Redfleld system and the 
open-systems Casida equations 

We now consider the open-systems LR-TDDFT formalism 
developed in section II, but applied to the Redfield master 
equation. The density-density response function of the Kohn- 
Sham-Redfield system is given by 



,[/?^(r),7(0] + ^7(f 



(55) 



We can now write Eq. |55]in a basis of Kohn-Sham-Mermin 
orbitals as 



1 



to + ^ks - 



(60) 

Using Eqs. 53-56 and inserting a complete set of Kohn- 
Sham-Mermin states, one obtaines the sum-over-states ex- 
pression. 



Jjijit) = -'E{'^'''(0«7<:;(0-:)^l(f)/j'''(f)-t;} + E'i;«%/(0 

k ijkl 

(56) 

We choose r^jj^,, to have the form of the Redfield tensor, but 
written in terms of Kohn-Sham-Mermin orbitals and eigenval- 
ues. 



(/|«(rO|;)(;|«(r)|/) | 



(/|«(r)|;)(;|«(rO|/) 



(61) 
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Eq.|5T|and Eq. |6T]are related through the Dyson-like rela- 
tion given in Eq. [36] To extract the poles of the interacting 
density-density response function in Eq. |5T| from that of the 
Kohn-Sham system in Eq. [61] a pseudo-eigenvalue equation 
must be solved for the squares of the complex transition fre- 
quencies. 



|(o2-n(co)|F = 0. 



(62) 



The operator Q.{co) can be written as a matrix in a basis of 
Kohn-Sham molecular orbitals (assuming a closed shell sys- 
tem) as 



^^{fi-fj){o^t + ^j)K.jkii(o) V^(A-//)«+A^fp) 

The explicit derivation of Eq.|62]and Eq.|63]is given in ap- 
pendix A. In Eq.|63] 

+fT'y''\-'^f'\{r,r'-co)\w)W) (64) 



IV. APPLICATION - SPECTRUM OF A C^+ ATOM FROM 
THE REDFIELD MASTER EQUATION 

As a simple demonstration, in this section we calculate the 
absorption spectrum of an atom in an optical resonator using 



the Redfield master equation 112 11 16111. The modes of the radia 



tion field act as a bosonic bath, leading to decay of the atomic 
excitations due to spontaneous and stimulated emission |24- 
3111 . We consider a hypothetical experimental setup similar 
to that used in 112 Sfl . The effect of the optical resonator is to 
modify the density of radiation field modes relative to the vac- 
uum, leading to an enhancement of the decay rate of atomic 
excitations at certain frequencies i28l roll l62ll . Although we 
do not have access to experimental data for in an optical 
resonator cavity, we use accurate experimental data taken in 
vacuum, together with a chosen cavity geometry to construct 
a numerically exact spectrum I16O1I . 



Spontaneous emission from tlie Redfield master equation 

For a single atom with zero center of mass velocity con- 
tained in an optical resonator cavity, the system-bath interac- 
tion is 



(68) 



and the bare Kohn-Sham linewidths and Lamb shifts are 
pven by the relation 



lA- 



(65) 



as in Eq |52l In principle, with the exact functional 
Jxc"'[n'^'' ,M,r''^], the exact poles of Eq.|5T|are recovered by 
solving Eq.|62l The operator Q.{co) is non-Hermitian, giving 
rise to complex eigenvalues corresponding to broadened exci- 
tation spectra. Q.{co) is also frequency-dependent and imagi- 
nary, both explicitly through the third term in Eq.|63]and im- 
plicitly through fxf™ in the coupling matrix Kijiii{co). The 
explicit frequency-dependence arises because the bare Kohn- 
Sham transitions are already broadened, even in the absence 
of Hartree-exchange-correlation effects. This is most easily 
seen by setting Kiji^i{o}) = in Eq.|63] Q.ijki is then diagonal 
and Eq. |62]reduces to a set of uncoupled equations given by 



Here, /i is the dipole operator for the atom, e,- and O); are 
respectively the polarization vector and frequency of the ith 
mode of the radiation field and V is the quantization volume. 
Qi and a] respectively destroy and create a photon in the ith 
mode of the cavity. The photon reservoir Hamiltonian is 



Hr ^Y^(£)ia\ai. 



(69) 



With the system-bath interaction and bath hamilto- 
nian specified, the Redfield tensor can be explicitly con- 
structed ieili . 

The entire atom-field system is taken to be in thermal equi- 
librium at inverse temperature j3, such that (Bqi ^ ^. With this 
condition, the atom can be assumed to be in it's groundstate 
and the effect of stimulated emission is neglected. We then 
need to only construct the matrix elements RaQaO appearing in 
Eq. |5T| For the real part of ^aOao one finds 



{« + A^)2 + (rfj)2-2,a)r^}^;^ = (o'Fij. (66) 
These are solved with the quadratic formula to yield 



aO — 



{2nf 



\JXao-ei?(Okg{(Ok,k)5[cOa(i - C0k)dQ.kd(0k. 

(70) 



The imaginary part of 7?onOa is given by 



c.^-,r-±(a)-+Af;), (67) A„o = ^L^ / j \%,rh\^^^sMd^,do.., 

y „ J J COan-COk 

which are precisely the poles of Eq. [61] (71) 
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where !^ denotes the principle value integral |'6l'l. is 
the frequency of a photon whose wave-vector magnitude is 
k— 1^1 . (S)a„ = Ea — E„ is the difference in atomic energy levels 
and /!„„ are matrix elements of the dipole operator between 
atomic wavefunctions. g{(£)) is the density of field modes in 
the cavity. In free-space, the density of field modes takes the 
formgf,-ee{(o) = fj^i,.i ■ This gives rise to a free-space natural 
linewidth 



-^free 
aO 



3c3 



(72) 



If one considers an experimental setup such as that used 
in 12811 . the density of field modes is modified to 



p{C0)=gfree{C0)M{C0), 

where the function M{co) is given by 



M{co) 



VT+F 



The linewidth in the cavity is then modified to 



aO — 



3c3 



M(a), 



>aO) 



(73) 



(74) 



(75) 



where we have neglected cavity edge effects and angular 
dependence of p. In Eq. |74] L is the length of the cavity and 



F = 



4R 



where R is the reflection coefficient of the cavity 



walls. By changing the mirror reflectivity and cavity length, 
suppression or enhancement of the spontaneous emission rate 
is possible. For our calculations, we choose cavity param- 
eters of R = 0.998 and L — 4.88 cm, leading to an overall 
enhancement. Using experimental data for the atomic energy 
levels and the transition dipole matrix elements of taken 
from [60], together with the specified parameters for the cav- 
ity geometry, we can explicitly construct Eq. |75] We include 
in our calculation the 3 lowest dipole allowed transitions of 
C^"*", which are li^2i^ — > ls^2s{2p,3p,4-p). The numerical 
values of the linewidths (imaginary part of the frequency) cal- 
culated in Eq. |75] are given in the fourth column of Table HI] 
Due to lack of experimental data on all atomic transitions, the 
Lamb shifts in Eq.|7T|cannot be explicitly evaluated. They are, 
however, estimated to be several orders of magnitude smaller 
and will be neglected in the following analysis. From the tran- 
sition dipole matrix elements we can also construct the oscilla- 
tor strengths. The "exact" response function constructed with 
these parameters is included in Figures 1-3. 



OQS-ATDDFT calculation of the spectrum of 

As a first step, we solve Eq.|62]using only an adiabatic ap- 
proximation to the exchange-correlation kernel. This corre- 
sponds to including the first term in Eq. |37] but entirely ne- 
glecting f^"''' (r , r' ; (B ) . 



For our calculations, we obtain the Kohn-Sham parameters 
of C^^ to be inputted in Eq.|63] using the real-space TDDFT 
package Octopus 063146511 . First, a ground-state DFT cal- 
culation is performed using the local density approximation 
(LDA) with the modified Perdew-Zunger (PZ) parameteriza- 
tion of the correlation energy ll66ll . For all calculations, the li^ 



core is replaced by a Troullier-Martins pseudopotential Il83ll . 
The Kohn-Sham eigenvalues and dipole matrix elements be- 
tween Kohn-Sham orbitals are computed, and substituted into 
Eq. 



I to obtain the bare Kohn-Sham linewidths, Pfj. These 



correspond to the real part of the Kohn-Sham-Redfield tensor 
(imaginary part of the bare Kohn-Sham frequency) of Eq. |57] 
and are given in the second column of Table |II] From the 
dipole matrix elements between Kohn-Sham orbitals, the bare 
Kohn-Sham oscillator strengths can be constructed. The real 
and imaginary parts of the bare Kohn-Sham response function 
constructed with these quantities are plotted in figures 1-3. 

Next, we perform a standard LR-ATDDFT calculation to 
obtain the matrix elements of the adiabatic Hartree-exchange- 
correlation kernel to be inputted in Eq. |63] The matrix ele- 
ments of the kernel obtained are given in table |III1 We also 
include the energies obtained from the standard LR-ATDDFT 
calculation in column 4 of Table H) For consistency, we 
use the LDA with modified PZ functional for the exchange- 
correlation kernel as well. 

Since the operator in Eq. |63] is explicitly frequency- 
dependent even when using an adiabatic kernel, Eq.|62] rep- 
resents a non-linear eigenvalue problem. We solye^it using the 
generalized eigenvalue algorithm presented in iItiI l72tl . The 
real part of the solutions to Eq. |62] are given in column 3 of 
Table U while the imaginary part is given in column 3 of Ta- 
ble In] The real and imaginary parts of the response function 
obtained are plotted in figures 1-3. 



Table I: Real part of the 3 lowest transition frequencies for in an 
optical resonator in a.u. 

Transition Bare Kohn-Sham OQS-ATDDFT ATDDFT Exact 



2s — !> 2p 
2s 3p 
2s 4p 



0.311 
1.116 
1.368 



0.443 
1.107 
1.361 



0.443 0.467 
1.107 1.180 
1.363 1.470 



Table II: Imaginary part of the 3 lowest transition frequencies for 
in an optical resonator in a.u. The last column includes the GL 
perturbation correction to the 2s 2p transition. 

Transition Bare Kohn-Sham OQS-ATDDFT Exact OQS-ATDDFT + GL 



2s 2p 
2s 3p 
2s 4p 



1.329 xlO 
5.162 xlO 
2.443 X 10- 



-2 
-2 



1.331 xlQ-^ 2.932 x 10 



5.159 xlO 



-2 



5.740 X 10 



-2 
-2 



1.805 XlO 



-2 



2.444 x10-3 1.089 x10-2 
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Figure 1 : Absorption Spectrum of including the 3 lowest dipole 
allowed transitions. The curves shown are: a) The bare Kohn- 
Sham spectrum (green-dashed), b) The spectrum obtained by solving 
Eq.|62]with an adiabatic exchange-correlation kernel (blue-solid), c) 
The numerically exact spectrum obtained using experimental data 
(red-dashed). Also shown is the stick spectrum obtained by solving 
the usual Casida equations for using ALDA (black-dotted). 



Figure 3: Real part of the density-density response function of C 
including the 3 lowest dipole allowed transitions. The effect of the 
photon bath is to broaden the dispersion over multiple frequencies. 
The curves shown are: a) The bare Kohn-Sham dispersion (green- 
dashed), b) The dispersion relation obtained by solving Eq|62|with 
an adiabatic exchange-correlation kernel (blue-solid), c) The disper- 
sion relation obtained using experimental data (red-dashed). 



-OQS-ATDDFT 
Bare Kohn-Sham absorption 
Exact 




(i)(a.U.) 



A 1.5 



Figure 2: Same as Figure 1, but with a close-up view of the 2s - 
and 2s 4p transitions. 



V. ANALYSIS 



Effect of using an adiabatic approximation to 



■3p 



Sham linewidths requires an additional frequency-dependent 
functional. This justifies a posteriori our separation of the 
exchange-coiTelation kernel in Eq. [37] into an adiabatic part 
and a frequency-dependent part due exclusively to bath ef- 
fects. 

To understand this situation better, we consider a "small 
matrix approximation" (SMA), in which a single occupied to 
unoccupied Kohn-Sham transition, i j, is completely iso- 
lated from all other transitions ll67l - l70ll . This is valid when 
the transition of interest is weakly coupled to all other excita- 
tions. In this case, Eq.l62lreduces to 



which is equivalent to the polynomial equation 



(76) 



(0^ + 2iTf](o- 



(77) 

If one assumes the adiabatic approximation, the coupling 
matrix is frequency independent, i.e. 



= 0. 



From Tables |T] and [III it is clear that using an adiabatic 
kernel in Eq. |62] gives essentially the same corrections to 
the real part of the energy as usual LR-ATDDFT, while giv- 
ing almost no correction to the imaginary part. This means 
that using an adiabatic kernel in OQS-TDDFT is expected to 
give the same reliable correction to the location of absorption 
peaks as usual LR-ATDDFT, while correcting the bare Kohn- 



K,jjj{(o)KK,j^ij 



(78) 



Eq. [TT] then reduces to a simple quadratic equation, with 
solutions given by 



CO = 



ir^^±Jico^J + /^^)^+4ico^ + A>;^)K.jjj. (79) 



11 



This shows that the eigenvalues retain their bare Kohn- 
Sham linewidths. It is evident from Table |III] that in gen- 
eral, the diagonal matrix elements of the kernel are apprecia- 
bly larger than the off diagonal elements, suggesting that the 
transitions are weakly coupled. A more rigorous criterion for 
the validity of the SMA is given in Eq. (11) of ||67|l- It can 
be shown that this criterion does in fact hold for the 3 lowest 
transitions of which we have included in our calculation. 



Ka)) = ^aV), 
(=1 

as well as the corresponding energies 



(81) 



(82) 



Table III: Matrix elements of the adiabatic Hartree-exchange- 
correlation kernel (Eq.|64) used in solving Eq.|62| 



Kernel matrix element Numerical value 



K2slp,2s2p 


8.035 xlO" 


2 




-4.827 X 10 


-3 


K2s4p,2s4p 


-2.528 X 10 


-3 


K2s2p.2s'ip 


1.031 xlO" 


2 


K2s2p.2s4p 


1.945 xlO" 


3 


J^2s3p,2s4p 


-4.730 xlO 


-4 



We also note that for the 2s 3p and 2s — !• 4p transitions, 
ATDDFT seems to provide much less reliable results then for 
the 2s ^ 2p transition. This is most likely related to the fact 
that we have truncated the occupied-unoccupied space to only 
the 3 lowest transitions. In the following analysis, we will 
focus on the low-lying 2s -^2p transition. 



Beyond an adiabatic approximation to /v/™ 

In the previous section, we discussed the effect of ap- 
proximating /"/'^"(r,r'; 0)) with an adiabatic (frequency in- 
dependent) kernel. We now investigate what the frequency- 
dependent contribution /^"'''(r,r';a)) must be to correct the 
bare Kohn-Sham linewidths. To formulate our construction 
rigorously, we examine the difference between the Kohn- 
Sham-Redfield tensor in Eq. [57] and the interacting Redfield 
tensor in Eq. |44] The bath correlation functions in both ex- 
pressions are the same. The difference lies in the eigenen- 
ergies and wavefunctions of the system Hamiltonain used to 
construct the Redfield tensor These are Kohn-Sham quan- 
tities in the first case and many-body quantities in the latter 
case. This suggests that the interacting Redfield tensor can 
be expanded in a Gorling-Levy (G-L) perturbation series in 
the electron-electron interaction coupling constant a, with the 
Kohn-Sham-Redfield tensor entering as the zeroth-order term 
in the series 1 73 - T^l : 



2„2 



Rabcdia) w Rahcd{Q) + aRabcd + <y- 



abed 



(80) 



For linear response from the ground-state as we consider 
here, we need only construct the quantities {^aOuo(o:)}. This 
is done by first expanding the ground-state and excited-state 
wavefunctions in a G-L perturbation series in a. 



These expansions are then substituted into the general ex- 
pression for the Redfield tensor to construct an expansion of 
the Redfield tensor at coupling constant a. For a = 1, we re- 
cover the interacting Redfield tensor in Eq. |44l since Eq. [81] 
and Eq. |82] then refer to wavefunctions and energies of the 
interacting system. For a = 0, we obtain a Redfield ten- 
sor /?aOao(0) written in terms of Kohn-Sham ground-state and 
excited Slater determinants, {|fl(0))}, and corresponding en- 
ergies {^^((O)}, which lie on the adiabatic connection with 
the interacting wavefunctions {|a(l))} and energies {^(((l)}. 
Due to the one-body nature of the system-bath coupling, only 
matrix elements of the tensor /?ao«o(0) containing singly- 
excited Kohn-Sham Slater determinants are non-zero. The 
zeroth-order term in the expansion of Eq.[80]then reduces to. 



RaQa0{Q) = R 



.ks 
•jij ■ 



(83) 



Here, the indices i, j label a pair of Kohn-Sham orbitals, in 
which an orbital 0,- occupied in the Kohn-Sham groundstate is 
replaced by an orbital occupied in the singly excited deter- 
minant |fl(0)). With the G-L expansion of the Redfield tensor 
well formulated, one can rigorously construct a G-L expan- 
sion of Jxc ™(r,r'; co) using the same general procedure out- 
lined in fl^. 

For the specific example of we consider here with the 
Lamb shifts neglected, Eq. [80] amounts to expanding r„o, in a 
G-L series as 



r„o(a) 



ar, 



(84) 



and determining the corresponding corrections to the bare 
Kohn-Sham linewidth. This is done explicitly in Appendix 
B for the ls'^2s^ ls^2s2p transition within the SMA. The 
first-order correction Fj^ 2s found to be 



^2p,2s = --^M{e2s - e2p)(e2.v - £2pf 

2 

3. 



d^r^siy)r^p{r) 
X [{2s2s\2s2s) - {2s2s\2p2p) ~ {2s2p\2s2p)\ , (85) 
where (^2s and are ground-state Kohn-Sham orbitals and 



{ij\kl) 



d r ; 



(86) 
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-OQS-ATDDFT 
--Exact 

-OQS-ATDDFT + GL 




Figure 4: Correction to the bare Kohn-Sham linewidth to first-order 
in GL perturbation theory. 



The frequency-dependent matrix element of the kernel 
f^f^{ry;(o) to first-order in G-L perturbation theory is 
found to be (see Appendix B): 



A' 



bath 

2s2p.2.s2p 



[CO) 



(® + 'r|;,2,)(ri,,.)- (87) 



To calculate the numerical value of Eq. [85] we evaluated 
the Coulomb integrals of Eq. [86] using the spectral-element 
code previously presented in |86], taking the orbitals from 
Octopus as input. The Octopus orbitals were obtained on a 
uniform mesh with a grid-spacing of 0.1 a.u. The spectral- 
element code projects the orbitals onto a tensorial Chebyshev 
basis for the efficient evaluation of Eq. [86] and we chose all 
parameters conservatively such that the main source of error 
is in the original representation of the orbitals. 

The imaginary part of the frequency of the 2s —>-2p transi- 
tion obtained by solving Eq. [62] with Eq.[87]included is given 
in column 5 of Table [II] In Figure 4, the imaginary part of 
the density-density response function with Eq.[87linculded is 
plotted in the vicinity of the 2s — > 2p transition. Because we 
have derived Eq. [87] within the SMA, there is no correction 
to the oscillator strength 1 67 , T^l and the entire effect seen in 
Figure 4 is due to a change in the linewidth. The contribu- 
tion from /'^"'''(r,r'; o) provides significant improvement to 
the linewidth, but is still far from the exact value, suggesting 
that higher-order corrections are important. 

The G-L expansion procedure outlined above is expected to 
provide large corrections to the bare Kohn-Sham linewidths in 
systems with strong static correlations. In this case, the Kohn- 
Sham single Slater determinants are a fundamentally poor de- 
scription of the interacting wavefunctions, and so the bath will 
interact with the Kohn-Sham system in a very different way 
than with the interacting system. By including just the first- 
order correction, one can strongly mix in excited determinants 
and provide a large correction to the linewidths. However, in 



the example we have considered, no static correlation is 
present which might explain why one needs to go to higher or- 
ders in the GL expansion to obtain the exact linewidth. Also, 
we have used the SMA to derive Eq. [85] and Eq. [87] so we 
have only allowed mixing in of the first excited Kohn-Sham 
determinant (see Appendix B). 



VI. CONCLUSION AND OUTLOOK 

We have formulated a general framework of LR-TDDFT 
for many-body open quantum systems, which in principle 
gives access to environmentally broadened spectra in a strictly 
ab initio manner. Our treatment is most applicable to micro- 
scopically derived master equations, in which the system and 
bath are both treated by starting from an underlying micro- 
scopic Hamiltonian. As an example, we analyzed the micro- 
scopically derived Redfield master equation for an atom inter- 
acting with a photon bath. In this case, the bath correlation 
functions were very well characterized since they depended 
only on the free radiation field density of states and the cavity 
geometry. In some cases, the bath correlation functions may 
be more complicated and need to be treated in an approxi- 
mate way. In the case of a phonon bath, one would need to 
be able to calculate the phonon density of modes as well as 
electron-phonon couplings. This could be done by obtaining 
the vibrational normal modes with a usual DFT calculation, 
and feeding the couplings and phonon density of states into an 
open-systems LR-TDDFT calculation for the broadened spec- 
trum. This would be particularly applicable to describing the 
absorption spectrum of impurity molecules imbedded in a lat- 
tice. Similarly, for the case of chromophores in a protein envi- 
ronment, one could simulate the protein using classical molec- 
ular dynamics to obtain the spectral density and then feed it 
into an open-systems LR-TDDFT calculation to compute the 
absorption spectrum of the chromophore. A similar treatment 
could be applied to molecules in liquid environments. These 
directions will be explored in our future work. 

The Redfield master equation is Markovian and as a re- 
sult, the pole structure of the OQS response functions are sim- 
ple. For Markovian master equations, there are only two pa- 
rameters characterizing the absorption spectrum; the location 
of the peaks and their width. As a result, the OQS Casida 
equations have relatively simple frequency-dependence. For 
non-Markovian master equations, the memory kernel is non- 
local in time, which coiTesponds to a frequency-dependent 
self-energy in the OQS response functions. The resulting 
lineshapes are asymmetric and may have many non-zero mo- 
ments. This is expected to give rise to much more compli- 
cated frequency dependence in the exchange-correlation ker- 
nel, since it must correct all of these moments. A simple step 
would be to investigate the first-order non-Markovian correc- 
tion derived from the cumulant expansion of the memory ker- 
nel [84]. This will be investigated in future work. 

This paper has focused on homogeneous broadening of the 
spectrum, which arises when the time-scale of the bath is 



13 



much faster than that of the electrons. This is the case of 
interest in OQS-TDDFT, since it implies that the bath can 
induce relaxation and dephasing of the electronic degrees of 
freedom. In the other limit of inhomogeneous broadening, the 
bath is static relative to the time-scale of the electrons. In this 
case, the external potential due to the nuclei is distributed in 
different configurations due to the local environment, but no 
relaxation and dephasing takes place. Inhomogeneous broad- 
ening can be well captured by performing usual closed LR- 
TDDFT calculations for different static nuclear configurations 
and then ensemble averaging. A realistic spectrum usually 
consists of both broadening mechanisms. This can be well 
captured by performing an open systems LR-TDDFT calcu- 
lation for each static nuclear configuration and then ensemble 
averaging over the different configurations afterwards. 

To perform an open-systems LR-TDDFT calculation us- 
ing modern electronic structure codes, the greatest challenge 



is probably the implementation of the algorithm in 117 iL 17211 



for solving the complex and non-linear eigenvalue prob- 
lem. This algorithm becomes expensive for a large occupied- 
unoccupied space of Kohn-Sham orbitals and a self consistent 
procedure becomes more efficient. Work towards implement- 
ing these capabilities in a numerical library is currently under 
investigation. 
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Xijki{(0) = UPiEa)] £ < ^ ' ^ p 

{a\aj^ai\b){b\a^jai\a) | 



CO-CO.b + lRlbab 



(90) 



In this basis, the linear response of the interacting 1 -particle 
reduced density matrix at frequency co is 



SPiji^) = Y.Xijkiio))dv,^,{co)ki, (91) 

kl 

and the linear density response is given by 



5«(r, co)^Y, Sp.jHUrMr). (92) 

ij 

In the Kohn-Sham-Redfield system, the same density re- 
sponse is produced from 



5n{r,co) = Y,5Yijico)(l)j{r)^i{r), 

ij 



(93) 



where 
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S7ij{co) = Y,x^Jkii^)Svksico)ki, 

kl 



(94) 



is the response of the Kohn-Sham-Redfield one-particle re- 
duced density matrix and 



APPENDIX A - DERIVATION OF THE CASIDA EQUATIONS 
FOR THE REDFIELD MASTER EQUATION 

In this appendix we derive Eq. |62]and Eq.|63]for the open- 
systems Casida equations. In keeping with the original deriva- 
tion by Casida |8], we introduce second quantized creation 
and annihilation operators, aj and a,-, for a one-particle or- 
bital basis. In what follows, this will be taken to be a basis 
of real molecular orbitals from a ground or equilibrium-state 
Kohn-Sham calculation. An arbitrary one-particle operator O 
is represented in this basis as 



fj-fi 



Xijki{(o) ^ 5ik5j, . 



0)-0)ij + ir^^j. 



Eq. (|94b can then be written as 



We now write 



(95) 



S7.j{co) = \ f (5v,,,(«),-, + 5v,(a))„ + 5vr"( 



(96) 



•J 



(88) 



The interacting density-density response function in Eq. 
(ISTT i is given by 



X„„(r,r';«) = £ (^),-(r)(/)/(r)x,;H(«)(^),(r')(/)/(r'), (89) 

ijki 

where 



5v„{co)ij + 5vT"{coh = l^K,ju{co)5Yki{co), (97) 

kl 



where 



Kijki{(o) 



^\/rfV0r(r)0;(r)|^ 



+/r"[«^1(r,r';co)U,(r')0,(r'). 



(98) 
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We can then re-write Eq.|96las 



fkH, 

E 

kl 



OikOjl 7 Kijui 

Jl-Jk 



(0 



5Ykli(0) = dvextiCO)ij. 

(99) 



We separate the Kohn-Sham-Redfield matrix into real and 
imaginary parts 



A-iT B 
B A + iY 







)-(: 









\ 5Vext*i(0) J ' 

or by applying the unitary transformation 



57(a)) 
8y {(o) 



(107) 



r^,,=r^ + iA^, (100) 



and separate particle-hole and hole-particle contributions as 



in 



V2 V -1 1 



(108) 



L 

i<i.ft>f, 



SikSji 



fl -fk 



fl^ fk 



£ Kijik{C0)dyik{C0)^dVex,{C0)ij 
klJk>fl 



-Kijkiia) 



A+B ir+coC 
iF+coC A-B 

-lSm{5vex,{0})) j 



(109) 



5Yki{(o) Without loss of generality we assume the applied perturba- 
tion to be real. Since the molecular orbitals are taken to be 
(-jQjjeal, the density response can be calculated from 5Re(57((o)) 
alone. From Eq. 11091 we obtain 



kl.A>fl 



SikSji- 



fk-fi 



fk-fi 



- E Kjiki{(o)5Ykii(o)^Sv,x,i(o)ji- 

kl,fk>fl 

We now define the following matrices: 



^Kjiikico) 



Aijkiico) = dikdj 



-^-K,jki{(0), 

Jk-JI 



^ijkl - 



Sik5j,r1fi 



fk - fl 

Bijkiico) = -Kijikico), 



and 



Cijki 



fk - fl 



8Yik{(o) 



(102) 



-1 



{A + B)-{iY+(oC)\A-B\ [lY+mC) 5Re(57(«)) 



(103) 



(104) 



(106) 



The matrices Bjjki and Cuki have the same form as in 



Ajjki is similar, but includes a contribution due to the Lamb 
shift and Fjjki is a new term. 

We can combine Eq. llOll and Eq. 1 1021 into a single matrix 
equation. 



5Re(5v,,t)(«). 



We introduce the matrices 



and 



s = -c[A~b] 'c, 



n((o) = -5"2 [a+b] s-^, 



(110) 



(111) 



(112) 



which have the same form as in the usual Casida equations. 
Eq. ll lOl can then be inverted to obtain 



(105) 3ie{5Y{co))=S-^{co^-Q,{co) + S-^T[A-B] ^Tr^ 



The poles of the density-density response function are ob- 
tained when the operator in brackets vanishes. Defining 



Ci{co)=D.{(o)-S-^T[A~B] 'r5"2 
- ico{S-^TC-^S^- +S^C-^TS-^), (114) 

this is equivalent to solving the pseudo-eigenvalue equation 
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|co^-n(co)|F = 0. 



(115) 



Returning to the basis of Kohn-Sham molecular orbitals, 
the matrix representation of Q.{(o) is 



APPENDIX B - FIRST-ORDER GORLING-LEVY 
PERTURBATION CORRECTION TO THE LINEWIDTH OF 
THE 2s 2p TRANSITION OF 

In this appendix, we derive the first-order correction to the 
bare Kohn-Sham linewidth in Eq.|85]as well as the frequency- 
dependent functional which gives rise to this correction. Our 
treatment closely parallels that used by Gorling in deriving 
the exact-exchange kernel in ll75ll . In what follows, the Lamb 
shifts are neglected. We also make the assumption that the two 
electrons in the Is^ core are frozen, and their effect on the two 
valence electrons is taken into account with an effective po- 
tential. This is in fact the case for our numerical calculations, 
since we replace the Is^ core by a pseudopotential. We then 
effectively have a two-electron singlet in which the interacting 
ground-state is 



VA2,2(r,r')^(r,r'|2.^(a = l)), 
and the first excited state is 



(117) 



5n{a,r,co) 



2{E2,{a)-E2p{a)) 



{E2s{a) - E2p{a)Y - (o) + iT2p.2.s{a)Y 
X {2s^{a)\h{r)\2s2p{a)){2s2p{a)\8v{a,(o)\2s^{ayini) 

For a~\, this expression refers to the density response of 
the interacting system, while for a = 0, it describes that of the 
Kohn-Sham system. In particular, dv{a — 1,G)) = 8\'ext{(i>) 
while 5v{a = 0,0)) = 5v\(a)). 5n(o:,a),r) = 8n{r,co) is in- 
variant with respect to the coupling constant a. We now ex- 
pand both sides of Eq. ll21l in a taylor series in a and equate 
coefficients of equal powers of a on both sides. At zeroth- 
order in a, we recover the Kohn-Sham response equation 



5n{r,co) = 



2(e2p — £2.5) 



ie2p -e2.r- (CO + iT% 2.^^ 



d\'{2s^{0)\pir)\2s2pi0)){2s2pi0)\p{r')\2s\0))dvMm) 



To evaluate the first-order terms, we need the expansions of 
the wavefunctions, energies and linewidths to first order in a. 
The expansions of the wavefunctions and energies are given 
by standard GL perturbation theory: 



\V2si)^ „ „ \¥2s2p{'J)) (123) 

fc2i fc2p 



|^t^^, (V^2.(0)|v..-v.-v.|v^2.2.(0)) |^^^^(^)^ (124) 

£2p — t2i 



^2*2 = (r2.v2(0)l<'« 



-v/,-v,-|v/2P(0)> 



(125) 



V/2.v2p(r,r')=(r,r'|2i2/7(a = l)). 



(118) 



We denote the respective energies of these two states by 
-^21-2 = '£'2^.2 (oc = 1) and £'2.v2p = £'2i-2p(a = 1). The corre- 
sponding Kohn-Sham ground-state is 



EL2p = {V2s2p{0)\Vee - Vh - | V^2.v2p(0)) • (126) 

Since the ground-state is a spin singlet, = ^ ^ and all 
quantities can be expUcitly evaluated. We find 



<I>2,2(r,r') = (r,r'|2/(a = 0) = (/)2.v(r)fc.(r'), (119) (r,r'|v4,2) = -j-^ -{2s2s\2s2p) 

^\£2s £2p) 

and Kohn-Sham first excited State is x [024 (r)02p(r') + fe (r')02p(r)] = (127) 



and 

4>2.2p(r,r') = (r,r'|252p(a = 0) = i=(fe.(r)02p(r') + fc(r')02p(r)). 



(120) r'lvfl - \ = {2s2s\2s2p) 



y-^^! (r,r |v^2.v2;/ = 

The respective energies are 2e2.v = E2^2((X = 0) and £2.5 + * 2V2(e2.s — £2^) 

£2p = £2.v2p(«-0). X [fe,(r)02,(r')] =0. (128) 

Our starting point is the linear density response at coupling 

constant a, in the subspace spanned by the Ls^2.s^ ls^2s2p i.e. the first-order corrections to the wavefunctions vanish 

transition, since (2s2s|2s2p) = 0. 
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The first-order G-L correction to the energy is 



= {2s2s\2s2s) - {2s2s\2p2p) - {2s2p\2s2p). (129) 

We now obtain the first-order correction to the bare Kohn- 
Sham linewidth. The linewidth for the 2s^ — > 2s2p transition 
at coupling constant a is 



To construct the kernel, we now equate coefficients of first- 
order in a in Eq. 11211 The result is 



X^^{(o.r,r')5v\(oy) + h"r{oi,vy)5v,{(oy) 



(135) 

Here, (ft),r') is the coefficient of the first-order GL ex- 
pansion of the potential, obtained from 



4 1, , 

r2p,2.v(0:) = -3(-) («2.2,2.v2p(«)) 
X A^(«2.s2,2.s2p(«))l('/2.v2(«)lMlr2.v2p(a)>|^ (BO) 

Expanding in a taylor series in a, at zeroth-order we re- 
cover the Kohn-Sham linewidth. 



X M(e2.v - e2p)|(r2.2(0)IAIr2.2p(0))|2. 
In terms of Kohn-Sham orbitals this is, 



(131) 



X M(£2,-e2p)| j £/^r02.(r)r02p(r)|' 
At first-order in a we find 



1p.2. = -^(^)^(e2.v - e2„)^M(e2, - £2p) 



(132) 



5v'(a,a),r) w 5vs{co,r) + a5v'(a),r), (136) 



and 



/n'""(«,r,r') 



2<22,2p 



4(e2.v - £2p){{£2s - e2p)<22,2p - '^ip 2J((® + l^%,2s))) 



2s2p 

{{£2.-£2pr~iC0 + ir%^2s)'y 



X [2^2.s{r)hp{r)^2s{r')(p2p{r')] 



(137) 



We now separate out the part of h"''™{(0,r,r') which con 
tains the correction to the Kohn-Sham linewidth T^p 2s' 



I bath 



(«,r,r') 



8(e2.s-e2p)(rL.2,)(a) + in 



2p.2s> 



{{£2s-£2pY-{CO + im2.Yf 

X [fe.(r)02p(r)fe.(r')02p(r')] ■ (138) 
The functional arising from this correction is then given by: 



X (r2.2(0)|Alr2.v2p(0)) (r2.v2(0)|AIV^2.2p) + (r2.2 I A I fe2p(0)) 



- 4(-)^ (£2. - e2p)^M(e2, - £2^X2,2,2^ 

X l(r2.s2(0)|Alr2.2p(0))|2 

Using Eqs. 11271 - [T29l we can now evaluate Eq. 11331 We 
find that. 



-1 , 



/r^«,r,r')= dh"dh"'x'' {(0,ry')hT"'{(0,r"y") 



(133) 



(139) 



Here, x is the inverse of the Kohn-Sham response func- 
tion 



^ 2p,2.s 



137 



M(£2,. -e2p)(e2.v-e2p) 
2 



/'(«,r,r') 



4(e2. - e2p) 



r02i(r)r02p(r) 
X [{2s2s\2s2s) - {2s2s\2p2p) - {2s2p\2s2p)] . (134) 

which is the correction to the bare Kohn-Sham linewidth we 
have included in Eq.[85] We can now ask what the frequency- 
dependent exchange-correlation kernel is which gives rise to 
this correction. This is in some sense a generalization of the 
exact-exchnage kernel to OQS, in that it is correct to first- 
order in G-L perturbation theory. However, the form of this 
functional will now depend on the bath. 



(e2.-e2p)2-(a) + ;r|;2j)2 

X [fe(r)(?)2p(r)fe.(r')02p(r')] (140) 

in the restricted space. For the open-systems Casida equa- 
tions, we need the matrix element of A'l'^jp 2i2p 

((0)ofEq.[l39l 

which is given by 



K2"2p.2s2p{0^) 



2(e2i — e2p 



(« + ;r|i.2,)(ri )-(i4i) 



This is the matrix element given in Eq. [87] The remaining 
part of Eq. 11381 we have not separated out changes only the 
location of the 2s^ 2s2p transition and not the width. In the 
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calculations of section V, we have replaced this contribution 
to the kernel with an adiabatic functional in solving Eq. |62l 

To understand Eq. ll41l better. we consider the SMA equa- 
tion in the 2s -^2p subspace: 

«^ + 2!r|^2p«-{(e2.-e2p)^ 
+ (r^^.2p)^ + 4(e2.-e2p)^2.2p,2.2p(»)}=0. (142) 
If we separate out the adiabatic and bath parts as 



^2i2p.2.s2p(«) - K2s2p,2s2p - TT" —T (« + l^2p^2s) {^2p.2s) 

^\t-2s t2pj 

(143) 

and substitute into Eq. ll42l we get: 

0)2 + 2Kr|^.2p + r\,,.2> - { (£2. - e2p)' 
+ {^2s2pf' +^^2s2p^2p.2s+ ^i^2.s - £2p)K2s2p .2s2p'^ 

= 0. (144) 
The solutions are 



» = -'(r^^.2p+r^p,2J 

± (£2. - e2p)2 - (r^p,2.)2 +4(62. - e2p)^2.2p.2.siM5) 

Since the term ij'2p2sf' small relative to the AT- 

DDFT shift, the effect of ^'252^ 2i2p('^) ^ simple correction 
to the bare Kohn-Sham linewidth by Fj.^ 2s- 
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